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A  computational  approach  to  homogeneous  nucleation  is  proposed,  based  on  Eulerian 
description  of  the  gas  phase  coupled  with  a  Lagrangian  approach  to  the  cluster  phase.  A 
continuum,  Euler  /  Navier-Stokes  solver  VAC  is  used  to  model  the  gas  transport,  and  a 
kinetic  particle  solver  is  developed  in  this  work  to  simulate  cluster  nucleation  and  growth. 
The  new  model  was  found  to  reproduce  well  the  known  theoretical  dimer  formation  equilib¬ 
rium  constants  for  two  gases  under  consideration,  argon  and  water.  Reasonable  agreement 
between  computed  and  available  experimental  data  was  found  in  terminal  cluster  size  dis¬ 
tributions  in  nozzle  water  expansions  in  a  wide  range  of  stagnation  pressures.  The  proposed 
approach  was  found  to  be  orders  of  magnitude  faster  than  a  comparable  approach  based 
on  the  DSMC  method. 


I.  Introduction 

Homogeneous  condensation  plays  an  important  role  in  many  atmospheric  and  technological  processes, 
and  understanding  of  its  physical  mechanisms  and  dependencies  is  critical  for  a  number  of  engineering 
applications.  One  of  such  applications,  pertaining  to  post  boost  vehicle  operations  at  very  high  altitudes,  is 
related  to  thruster  plume  expansion  into  surrounding  rarefied  atmosphere.^  It  is  well  known  that  particulates 
of  different  kind  are  the  main  contributor  to  sunlight  scattering  observed  in  high  altitude  plumes.  The  effect 
of  sunlight  scattering  in  plumes  in  which  neither  carbon  soot  nor  alumina  particles  were  present,  with  the 
specific  example  of  the  Apollo  8  translunar  injection  burn,^  indicates  that  particles  must  be  formed  in  the 
rapid  expansion  of  the  exhaust  to  rarefied  atmosphere,  mostly  from  the  condensation  of  water  vapor  and 
other  combustion  products  in  the  plume. 

Condensation  in  the  rapidly  expanding  flows  has  been  observed  experimentally  as  early  as  mid-30s,^  and 
has  been  extensively  studied  in  the  following  decades  (see  for  example  Ref.  4  and  the  references  therein). 
Computational  modeling  of  expanding  condensing  flows  has  a  shorter,  although  still  a  respectable  history. 
In  the  past,  two  different  approaches  have  been  used  to  describe  homogeneous  condensation  and,  in  par¬ 
ticular,  cluster  nucleation  (formation  of  small  clusters  from  monomers)  in  non-equilibrium  environment  of 
rapid  expansions.  In  the  first  approach,  based  on  the  classical  nucleation  theory  (CNT)®  and  equilibrium 
thermodynamics,  the  key  process  is  the  formation  of  the  smallest  stable  droplets  possible,  so-called  critical 
clusters,  through  unimolecular  reactions  of  cluster  growth  and  decay.  The  classical  theory  calculates  the 
condensation  and  evaporation  rates  using  the  Gibbs  distributions  and  the  principle  of  detailed  balance,  and 
the  nucleation  rate  is  then  calculated  assuming  a  steady  state  condition.® 

The  main  principles  of  the  classical  nucleation  theory  in  combination  with  the  conventional  compressible 
Navier-Stokes  gas  dynamic  equations  were  used  by  a  number  of  researchers  to  numerically  predict  multi- 
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dimensional  condensing  flows  (see,  for  example,  Refs.  7-9).  The  important  part  of  these  models  is  the 
creation  of  cluster  nuclei  at  some  critical  size  that  depends  on  local  gas  conditions.  The  nucleation  rate  is 
governed  by  CNT,  and  droplet  growth  can  be  derived  on  the  basis  of  heat  transfer  conditions  surrounding 
the  droplet  (the  description  of  Ref.  10  was  used  in  Ref.  8).  Both  Lagrangian®  and  Eulerian^  description  of 
condensed  droplets  was  used  in  the  literature. 

An  alternative  approach  to  modeling  homogeneous  condensation,  is  based  on  some  assumed  shape  of 
the  droplet  size  distribution  function,  usually  lognormal.  In  Ref.  11,  this  assumption  is  coupled  with  the  a 
modified  form  of  the  Hertz-Knudsen  equation,  which  gives  a  droplet-gas  mass  transfer  rate  as  the  difference 
between  incoming  fluxes  from  the  gas  phase  and  evaporative  fluxes  from  the  droplet;  a  standard  Eulerian 
description  was  used  to  model  the  two-phase  flow.  In  Ref.  12,  viscous  compressible  reduced  Navier-Stokes 
equations^®  are  used  for  the  gas  phase,  while  polydisperse  particle  behavior  is  described  by  an  Eulerian 
aerosol  moment  model  which  accounts  for  particle  transport  due  to  convection,  diffusion,  inertia  and  ther¬ 
mophoresis,  as  well  as  particle  dynamics  due  to  coagulation,  nucleation  and  condensation.  Yet  another 
numerical  approach,  which  uses  many  of  the  CNT  assumptions,  and  has  been  applied  mostly  to  turbulent 
condensing  flows,  is  based  on  a  semi-Lagrangian  treatment  of  droplets. Semi-Lagrangian  methods  com¬ 
bine  both  Eulerian  and  Lagrangian  points  of  view:  a  scalar  field  is  discretized  on  a  Eulerian  grid,  but  is 
advanced  in  time  using  a  Lagrangian  technique. 

While  different  methods  were  applied  to  predict  cluster  nucleation  and  growth  in  gas  flows,  most  of  the 
researchers  that  used  the  classical  nucleation  theory  applied  an  Eulerian  approach  to  the  gaseous  phase, 
usually  based  on  the  solution  of  full  or  reduced  Navier-Stokes  equations.  A  different  strategy  was  proposed 
in  Ref.  16,  where  a  particle-based  direct  simulation  Monte  Carlo  (DSMC)  method^^  was  used  to  compute 
the  gas  flow.  A  Lagrangian  technique  was  applied  to  model  cluster  evolution.  Similar  to  Ref.  7,  new  clusters 
were  created  at  a  critical  size,  and  their  further  growth  was  calculated  with  the  CNT  approach. 

Some  important  assumptions  of  the  CNT,  such  as  unimolecular  reactions  of  cluster  growth  and  decay  and 
the  use  of  the  principle  of  detailed  balance  that  implies  thermodynamic  equilibrium,  limit  its  applicability 
as  a  prediction  tool  for  highly  non-equilibrium  flows,  such  as  rapidly  expanding  plumes.  In  such  flows,  the 
impact  of  thermal  non-equilibrium  between  gas  and  particles  is  expected  to  significantly  impact  the  growth 
rates  and  cluster  size  distributions.  Moreover,  the  cluster  size  distribution  may  have  a  significantly  more 
complex  shape  than  the  lognormal  distribution  often  used  in  the  literature.  An  illustrative  example  of  such 
complex  distributions  was  provided  in  Ref.  18,  where  the  terminal  cluster  size  distributions  were  measured 
in  water  and  ammonia  expansions  for  a  wide  range  of  stagnation  pressures  and  temperatures;  the  results 
were  obtained  by  doping  the  water  and  ammonia  clusters  by  one  Na  atom,  which  was  photoionized  close  to 
the  threshold  without  fragmentation. 

The  experimental  study^®  showed  that  for  lower  pressures,  the  size  distribution  is  exponential;  for  higher 
pressures,  the  size  distribution  approaches  the  lognormal  profile,  and  for  intermediate  pressures,  it  was  a 
complex  bi-modal  shape.  The  transition  from  the  exponential  to  bi-modal  shape  was  explained  by  changing 
governing  mechanisms  of  cluster  growth.  For  lower  pressures,  the  clusters  grow  mostly  through  monomer 
sticking,  while  at  higher  pressures,  the  main  mechanism  is  cluster-cluster  collisions  and  coalescence.  The 
bimodal  shape  of  the  cluster  size  distribution  function  for  intermediate  plenum  pressures  was  attributed  in 
Ref.  18  to  the  processes  of  coalescence  of  small  particles  (such  as  dimers  and  trimers)  on  larger  clusters 
and  of  coagulation  of  larger  clusters.  A  bimodal  distribution  of  cluster  sizes  was  measured^®  for  a  number 
of  chamber  pressures,  varied  by  up  to  an  order  of  magnitude;  typically,  it  was  observed  when  the  average 
cluster  size  was  from  below  100  to  about  1000.  These  cluster  sizes  are  believed  to  be  largely  occurring  in  a 
number  of  applications,  including  rocket  thruster  plumes. 

The  inability  of  CNT  based  methods  to  accurately  predict  the  cluster  size  distributions  in  strongly 
nonequilibrium  flows  dictates  the  use  of  the  second  approach,  known  as  the  kinetic  approach,  which  treats 
nucleation  as  the  process  of  kinetic  chemical  aggregation.^®  Unlike  CNT,  the  kinetic  approach  does  not 
assume  local  thermodynamic  equilibrium.  Instead,  a  microscopic  process  of  the  interactions  of  monomers 
and  clusters  is  described  either  analytically  via  a  mathematical  model,  e.g.,  by  the  Smoluchowski  equations 
where  the  interaction  between  particles  is  modeled  by  the  reaction  rates,®®’ or  in  computer  simulations, 
e.g.,  in  molecular  dynamics  calculations  where  the  interaction  is  modeled  by  an  interaction  potential.®®’®®  It 
is  well  known  that  the  application  of  either  the  Smoluchowski  equations  or  the  molecular  dynamics  approach 
to  the  modeling  of  cluster  evolution  in  multi-dimensional  thruster  plume  flows  is  computationally  unfeasible. 

A  more  promising  direction  in  modeling  rapidly  expanding  condensing  flows  is  the  use  of  the  DSMC 
method.  As  a  numerical  approach  to  the  Boltzmann  equation,  it  is  applicable  to  a  large  range  of  flow  condi- 
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tions.  In  this  method,  cluster-cluster  and  cluster-monomer  interactions  including  the  multi-body  reactions 
of  cluster  nucleation  can  be  seamlessly  incorporated.  Over  the  last  several  years,  the  DSMC  method  has 
been  extensively  and  successfully  applied  to  modeling  the  processes  of  cluster  formation  and  evolution  in 
supersonic  jets.^^’^^  The  work  of  these  authors^®  extended  the  kinetic  dimer  formation  approach  of  Ref.  27, 
who  assumed  that  a  ternary  collision  always  results  in  a  dimer  formation,  to  include  molecular  dynamic 
(MD)  simulations  for  obtaining  information  on  the  probability  of  dimer  formation  in  such  ternary  collisions. 
The  work^®  used  a  temperature-dependent  probability  of  formation  of  argon  dimers.  Another  DSMC-based 
model,  which  treats  both  cluster  nucleation  and  evaporation  (RRK^®  technique  was  used  for  the  latter)  from 
the  principles  of  the  kinetic  theory,  was  introduced  in  Ref.  30. 

The  downside  of  using  the  DSMC  method  as  the  modeling  approach  for  condensing  plumes  is  its  high 
computational  cost.  It  may  be  applied  to  relatively  low  density  plumes,  when  the  typical  size  of  clusters  does 
not  exceed  lOO-mers.  For  higher  pressures,  this  approach  becomes  prohibitively  expensive.  The  most  serious 
numerical  limitation  of  the  DSMC  method  is  related  to  the  fact  that  a  large  number  of  simulated  particles 
has  to  be  computed.  The  required  number  of  simulated  particles  generally  increases  as  for  2-D  problems, 
and  n®  for  3-D  problems,  where  n  is  the  gas  density.  Most  of  the  simulated  particles  are  monomers;  the 
statistical  scatter  for  cluster  species  is  therefore  extremely  high  as  compared  to  the  monomers.  The  use  of 
species  weights  in  the  DSMC  method  is  questionable,  since  clusters,  especially  for  larger  pressures,  are  not 
a  trace  species,  and  thus  strongly  impact  the  flow  properties  through  the  heat  release  during  the  nucleation 
and  cluster  growth  process. 

The  main  objective  of  this  work  is  the  development  of  a  new  method,  that  would  combine  the  com¬ 
putational  efficiency  of  an  Eulerian  continuum  approach  and  the  physical  accuracy  of  a  Lagrangian  kinetic 
approach.  The  proposed  method  integrates  an  Eulerian  approach  for  monomer  gas  flow  based  on  the  solution 
of  Navier-Stokes  equations,  with  a  Lagrangian  approach  for  clusters  based  on  a  DSMC-like  particle-based 
algorithm.  The  work  is  built  on  the  previous  effort®®  where  the  first-principles  model  of  homogeneous  con¬ 
densation  was  formulated,  and  all  of  the  most  important  processes  of  cluster  nucleation  and  evolution  were 
considered  at  the  microscopic  level.  The  processes  included  in  the  model®®  are  (i)  creation  of  dimers  through 
the  collision  stabilization  of  collision  complexes,  (ii)  elastic  monomer-cluster  collisions  that  change  the  trans¬ 
lational  and  internal  energies  of  colliding  particles,  (iii)  inelastic  monomer-cluster  collisions  that  result  in 
monomer  sticking,  (iv)  cluster-cluster  coalescence,  (iv)  evaporation  of  monomers  from  clusters.  All  these 
processes  are  present  in  the  new  method.  In  the  next  sections,  the  details  of  the  method  are  discussed,  and 
the  homogeneous  nucleation  rates  in  argon  and  water  thermal  bath  environments  are  analyzed,  followed  by 
the  validation  study  that  focuses  on  comparison  of  cluster  growth  in  plumes  with  available  experimental 
data  on  terminal  cluster  size  distribution. 


II.  Numerical  approach 

The  main  idea  of  the  present  numerical  method  is  to  calculate  gas  flow  solving  the  compressible  Navier- 
Stokes  equation,  model  the  nucleation  process  starting  from  the  dimer  formation  and  up  using  the  elementary 
kinetic  theory  for  cluster-cluster  and  cluster-monomer  collisions,  and  and  exchange  the  information  between 
the  continuum  and  kinetic  parts  of  the  simulation  through  source  terms,  so  that  these  parts  are  fully  coupled. 
Similar  to  the  DSMC  method,  a  finite  number  of  simulated  clusters  replace  the  real  ones,  so  that  each 
simulated  cluster  represents  a  large  number  of  real  particles. 


A.  Eulerian  approach  to  gas  phase 


The  Eulerian-Lagrangian  approach  with  a  two-way  coupling  developed  in  Ref.  31  to  model  two-phase  plume 
flows  represents  the  computational  framework  of  the  new  condensation  model.  Gas  properties  are  computed 
using  an  Eulerian  approach  based  on  the  solution  of  the  Navier-Stokes  equations  with  appropriate  source 
terms  that  take  into  account  the  impact  of  condensation  process  and  clusters  on  the  gas  flow, 


P 


DY 

~LH 


Dp 

Dt 


+  P(VV)  =  M, 


DY  ^  ^ 

p~j^  V  ■  —  Di, 

V  +  VP  ■  V  -  V  Tij)  ■  V 


Q, 


3  of  14 

Distribution  A:  Approved  for  public  release;  distribution  unlimited 


where  M,  Di,  and  Q  are  the  corresponding  mass,  momentum,  and  energy  source  terms  that  define  the  impact 
of  condensation  on  gas  molecules. 

The  Navier-Stokes  equations  are  solved  using  Versatile  Advection  Code  (VAC)^^  modified  to  include 
the  above  particle  source  terms.  Particle  properties  are  determined  by  Lagrangian  tracking  of  particles 
through  the  gas  flowfield  and  statistical  averaging  of  particle  parameters.  For  the  gas  phase,  an  explicit  time 
integration  is  used,  and  the  Navier-Stokes  equations  are  solved  using  the  TVD-Lax-Friedrichs  scheme  with 
minmod  limiter.  For  the  particle  phase,  a  fourth  order  Adams-Moulton  method  is  used  to  integrate  particle 
equations  of  motion. 

In  the  current  implementation,  the  clusters  are  assumed  to  be  in  the  translational  equilibrium  with  gas, 
that  is,  their  macroscopic  velocity  and  translational  temperature  are  assumed  to  be  equal  to  the  corresponding 
parameters  of  the  gas.  At  each  time  step  At,  the  clusters  are  moved  by  v^At,  where  Vj  is  the  velocity  of  the 
t-th  cluster.  Then,  the  cluster  collision  relaxation  processes  are  modeled  at  the  kinetic  level.  These  processes, 
that  include  the  formation  of  new  dimers,  monomer-cluster  collisions  that  involve  energy  transfer  between 
internal  and  translational  modes  of  colliders,  cluster-cluster  coalescence,  and  cluster  growth  and  shrinking  due 
to  monomer  sticking  and  evaporation,  are  described  in  detail  below.  After  the  cluster  relaxation  processes, 
the  changes  in  cluster  mass  and  internal  energy  is  evaluated,  and  then  used  to  calculate  the  RHS  of  the 
Navier-Stokes  equations. 

B.  Lagrangian  approach  to  cluster  formation  and  evolution 

Replacing  the  kinetic  modeling  of  gas  transport  with  a  continuum  approach  is  justified  by  the  proximity  to 
equilibrium  of  velocity  distribution  functions  of  gas  molecules  in  condensing  plumes,  where  the  gas  density  is 
fairly  large,  and  mean  free  path  is  many  orders  of  magnitude  smaller  than  the  characteristic  flow  size  (usually 
nozzle  throat  or  exit  diameter).  The  cluster  nucleation  and  evaporation  processes,  though,  require  kinetic 
treatment  for  a  number  of  reasons,  most  notably  non-equilibrium  cluster  size  distribution  and  the  departure 
from  equilibrium  of  cluster  internal  energies.  Such  a  kinetic,  Lagrangian  treatment  is  therefore  proposed  in 
the  present  work.  Although  the  Lagrangian  approach  to  modeling  cluster  nucleation  follows  to  some  extent 
the  first-principle,  fully  kinetic  approach  of  the  previous  work,^°  it  has  a  number  of  key  differences,  mostly 
related  to  the  fact  that  monomers  are  simulated  at  the  continuum  level,  and  thus  some  approximation  has 
to  be  used  to  include  cluster-related  collisions  that  involve  monomers. 


1.  Dimer  formation 


Dimers  are  formed  as  a  result  of  a  collisional  stabilization  of  collision  complexes  consisting  of  two  monomers 
that  collide  with  third  particles  during  their  life  time;  the  third  particle  is  needed  to  carry  away  extra  energy 
and  thus  stabilize  the  dimer.  In  each  stabilization  event,  there  is  also  an  energy  release  from  the  potential 
levels  of  two  monomers  that  formed  the  collision  complex,  to  the  internal  energy  modes  of  the  newly  created 
dimer  and  the  third  particle,  and  the  translational  modes  of  their  relative  motion. 

The  change  in  the  dimer  number  density  as  a  result  of  new  dimers  formed  in  each  cell  over  a  single  time 
step  And  is  calculated  from  the  known  recombination  rate  and  macroscopic  gas  properties  in  the  cell  as 
follows 

And  =  Krecn^At.  (1) 


Here,  Krec  is  the  recombination  rate  and  n  is  gas  number  density.  Then,  the  number  of  newly  created 
dimers  is  given  by  where  Fnum  is  the  number  of  real  clusters  represented  by  one  simulated  cluster 

(similar  to  Fnum  traditionally  used  in  the  DSMC  method)  and  14  is  the  cell  volume.  Generally,  any  form  of 
temperature  dependence  may  be  used  to  define  Krec]  in  this  work,  a  temperature  dependence  similar  to  the 
well  known  Arrhenius  dependence  is  used,  KreciT)  =  A  x  exp(— CT).  In  this  equation,  constants  A,  B 
and  C  may  be  chosen  either  from  values  known  in  the  literature,  or  selected  to  reproduce  analytical  dimer 
formation  equilibrium  constants. 

The  initial  position  of  each  formed  dimer  is  selected  uniformly  within  the  cell,  its  initial  velocity  is  set 
equal  to  the  macroscopic  velocity  of  the  gas  in  the  cell,  and  the  initial  cluster  internal  energy  is  sampled 
as  follows.  First,  the  total  available  energy  in  the  collision  complex-third  particle  collision  is  assumed  to  be 
equal  to 


^tot  —  (  2 


4  -  2ar 


4 -2a, 


kT, 


(2) 
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where  ^int  is  the  number  of  internal  degrees  of  freedom  of  the  monomer  (zero  if  atom),  and  a  is  the  VHS 
model  parameter.  In  this  expression,  the  4  —  2am  term  corresponds  to  the  number  of  relative  translational 
degrees  of  freedom  in  the  monomer-monomer  collision,  and  am  corresponds  to  this  type  of  collision.  Similarly, 
4  —  2ac  is  the  number  of  relative  translational  degrees  of  freedom  in  the  collision  complex-monomer  collision. 
After  that,  the  total  energy  is  increased  by  evaporation  (dimer  dissociation)  energy  Eevap,  and  then  split 
between  the  newly  created  dimer  and  the  third  particle  using  the  Larsen-Borgnakke  (LB)^^  procedure. 
This  procedure,  initially  developed  to  model  energy  transfer  between  translational  and  rotational  modes  of 
colliding  molecules,  is  based  on  the  assumption  that  after-collision  relative  translational  and  internal  energy 
modes  will  be  populated  according  to  the  local  equilibrium  distribution  functions.  In  this  work,  the  energy 
transfer  was  assumed  to  include  all  available  energy  modes,  i.e.  energy  of  relative  motion  of  the  dimer- 
monomer  pair,  the  internal  energy  of  the  monomer,  and  the  internal  energy  of  the  dimer.  For  monomers, 
only  rotational  modes  are  assumed  to  be  exited,  since  at  low  gas  temperatures  in  expanding  plumes  (below 
300  K  in  this  work)  the  excitation  of  vibrational  modes  is  negligible.  The  number  of  internal  degrees  of  the 
dimer  is  calculated  from  the  dimer  heat  capacity  Cy 

-  3,  (3) 

where  n  is  the  number  of  monomers  in  the  clusters  (n  =  2  for  dimer).  Note  that  this  expression  is  also 
utilized  for  larger  clusters.  The  values  of  heat  capacities  used  in  this  work  for  argon  and  water  clusters  are 
taken  from  Ref.  30. 

The  dimer  formation  procedure  thus  results  in  the  formation  of  dimers  at  a  given  temperature  dependent 
rate,  and  each  of  these  dimers  is  characterized  by  unique  internal  energy  that  is  then  used  in  the  cluster 
collision  and  evaporation  processes. 

2.  Inelastic  monomer- cluster  collisions 

The  interaction  between  monomers  and  clusters  is  an  important  process  that  results  in  energy  transfer 
between  the  internal  energy  modes  of  clusters  and  translational  modes  of  colliders.  The  change  in  the 
cluster  internal  energies  greatly  impacts  the  evaporation  rates,  and  thus  needs  to  be  properly  modeled.  It 
was  pointed  out  in  Ref.  30  that  when  the  LB  model  is  used  to  simulate  the  energy  transfer  in  monomer  - 
cluster  collisions,  it  is  reasonable  to  introduce  an  inelastic  collision  relaxation  number  Z,  which  defines  the 
probability  that  a  cluster  will  experience  an  inelastic  collision  leading  to  a  change  in  its  internal  energy  in  a 
single  collision  as  Pi  =  ^.  This  means  that  only  every  one  out  of  every  Z  collisions  of  a  cluster  will  lead  to 
the  change  in  its  internal  energy.  In  every  collision  that  involves  such  a  change,  the  after-collision  energies 
are  selected  according  to  the  local  equilibrium  distribution  functions.  It  is  thus  similar  to  the  rotational 
and  vibrational  relaxation  numbers  Zy  and  Zy  widely  used  in  the  DSMC  method.  The  values  of  Z  were 
chose^°  to  provide  good  agreement  with  known  theoretical  dimer  formation  equilibrium  constants  for  argon 
and  water. 

A  similar  approach  is  used  in  the  present  model,  with  one  signihcant  exception.  Among  many  monomer- 
cluster  collisions,  only  those  that  cause  the  cluster  internal  energy  change  are  important  in  terms  of  cluster 
evolution,  and  thus  only  those  collisions  were  modeled.  Since  the  value  of  Z  is  typically  larger  than  one 
(for  example,  for  water  it  is  10  for  temperatures  between  100  K  and  300  K),  such  an  approach  allows  for 
significant  reduction  in  computational  time.  The  approach  therefore  is  reduced  to  the  analysis  of  each  cluster 
in  the  computational  domain  in  terms  of  possible  inelastic  collisions  with  monomers  as  follows. 

First,  note  that  the  probability  that  a  cluster  will  experience  an  inelastic  collision  leading  to  a  change 
in  its  internal  energy  during  time  r  is  equal  to  p  =  1  —  exp{vT  j Z) ,  where  n  is  collision  frequency.  Time  to 
the  next  inelastic  collision  can  be  sampled  as  Tie  =  ^  log(3i)^/j^,  where  3?  is  a  random  number  uniformly 
distributed  between  0  and  1.  The  algorithm  to  model  inelastic  collisions  for  a  given  cluster  over  time  step 
At  is  thus 

1.  Set  tiocai  to  0 

2.  Calculate  n  and  Z  as  functions  of  gas  macroparameters 

3.  Change  Uocai  =  kocai  -  \og{^)Z/v 

4.  If  tiocai  >  At,  go  to  the  next  cluster 


5  of  14 

Distribution  A:  Approved  for  public  release;  distribution  unlimited 


5.  Perform  inelastic  collision. 

6.  Go  to  step  3 

Unlike  Ref.  30,  where  the  kinetic,  microscopic  information  on  monomers  that  includes  their  individual 
energy  states  and  velocities  is  available,  the  present  method  provides  only  the  macroscopic  information  such  as 
temperature  and  number  density.  While  sufficient  to  calculate  the  local  collision  frequency  and  temperature- 
dependent  internal  energy  relaxation  number  Z,  this  is  not  enough  to  simulate  monomer-cluster  collisions 
(Step  5)  at  the  kinetic  level.  In  order  to  do  that,  an  additional  information  about  the  velocity  and  energy 
distribution  functions  of  monomers  is  necessary.  In  this  work,  the  internal  energy  of  the  colliding  monomer 
and  the  relative  translational  energy  of  the  colliding  monomer-cluster  are  sampled  from  the  corresponding 
equilibrium  distributions.  The  total  collisional  energy,  which  is  a  sum  of  these  two  energies  and  internal 
energy  of  the  cluster,  is  then  redistributed  between  the  relative  translational  and  the  internal  modes  of  the 
cluster  and  the  monomer  using  the  LB  model.  The  numbers  of  the  corresponding  degrees  of  freedom  are 
defined  as  described  in  the  dimer  formation  section. 


3.  Cluster  growth  and  evaporation 

The  key  processes  that  determine  small  cluster  evolution  is  sticking  and  evaporation  of  monomers  off  the 
clusters.  In  order  to  drastically  reduce  the  requirements  to  the  minimum  time  step  used  in  the  simulation,  and 
provide  accurate  account  of  evaporation  and  sticking  events  of  a  single  cluster,  the  growth  and  evaporation 
processes  are  combined  in  a  single  algorithm  as  follows. 

The  cluster  sticking  rate  is  calculated  as  Vs  =  nPs{acg),  where  Ps  is  the  probability  that  a  monomer 
will  stick  to  the  cluster  after  the  collision,  CTc  is  monomer-cluster  collision  cross-section  calculated  using  the 
hard  sphere  model,  and  g  is  relative  collision  velocity.  In  the  hard  sphere  model,  where  the  collision  cross 
section  is  written  as  Trd^,  the  collision  diameter  d  is  written  as  the  average  of  the  diameter  of  the  colliding 
monomer  taken  from  the  VHS  modeb^  and  the  cluster  diameter  obtained  through  an  empirical  correlation 
used  extensively  in  the  past  (see,  for  example.  Ref.  16), 

dc  =  2-{A-i^  +B),  (4) 


where  A  and  B  are  species-dependent  constants,  and  i  is  the  number  of  monomers  in  the  cluster.  In  this 
work,  the  values  of  A  and  B  were  2.3  x  m  and  3.4  x  10“^°  m  for  argon, and  1.9  x  10~^°  m  and 

2.4  X  10“^°  m  for  water. 

For  monomer-cluster  sticking  collision  probability,  an  empirical  dependence  of  the  probability  on  the 
species  diameter  d  and  mass  m  given  in  Ref.^®  is  used,  that  may  be  written  as 

"7  ) '  ,  (5) 

niN  +  mi  J 


d% 


(djv  +  di)"^ 


where  indices  N  and  1  refer  to  the  cluster  of  a  size  N  and  monomer,  respectively. 

To  evaluate  the  rate  of  evaporation  of  monomers  off  the  cluster  surface,  the  RRK  modeb®  is  used,  similar 
to  Ref.  30.  Following  Ref.,^®  the  evaporation  rate  ke  is  calculated  as 


ke  =  vNs 


-E 

h^int 


evap 


3N-7 


(6) 


Here,  N  is  the  number  of  monomers  in  the  cluster,  v  is  the  vibration  frequency,  Ng  is  the  number  of  surface 
atoms,  Eevap  is  the  evaporation  energy,  and  Eint  is  the  cluster  internal  energy.  For  dimers,  the  exponent 
3N  —  7  is  replaced  with  1.  The  number  of  surface  atoms  Ng  is  N  for  N  <  5,  N  —  1  for  4  <  N  <7,  and 
(367r)^/^(A^^/^  —  1)^  for  N  >  6.  The  vibration  frequency  was  taken  to  be  2.68  x  10^^  s“^  for  water  clusters, 
and  10^^  s“^  for  argon  clusters.^® 

With  the  evaporation  and  sticking  rates  defined  by  the  above  expressions,  the  following  algorithm  is  used 
to  model  sticking  and  evaporation  processes. 


1.  Set  tiocai  to  0. 

2.  Calculate  sticking  and  evaporation  rates  Vg  and  nue 
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3.  Change  tiocal  =  Uocal  -  log(5ft)/(l^e  +  l^s) 

4.  If  tiocal  >  proceed  to  the  next  cluster 

5.  If  3?  >  Vsli^s  +  ^'e))  increase  cluster  size  by  one  monomer  (see  below) 

6.  Otherwise,  evaporate  one  monomer  (see  below).  If,  after  evaporation,  the  cluster  becomes  a  monomer, 
remove  this  cluster  from  the  simulation  and  proceed  to  the  next  cluster. 

7.  Return  to  step  2 

For  cluster  growth,  the  monomer  internal  energy  and  relative  translational  energy  are  sampled  from  the 
corresponding  equilibrium  distributions,  and  the  after-sticking  cluster  internal  energy  is  equal  to  the  sum 
of  cluster  pre-collisional  internal  energy,  internal  energy  of  the  monomer,  relative  translational  energy,  and 
evaporation  energy  Ecvap- 

For  cluster  evaporation,  the  cluster  internal  energy  is  decreased  by  Ecvap,  and  the  remaining  energy  is 
redistributed  between  the  cluster  internal  modes,  internal  modes  of  the  departing  monomer,  and  relative 
translational  modes  using  the  LB  model.  Note  here  that  only  cluster  internal  energy  is  calculated,  while  the 
cluster  velocity  and  monomer  properties  are  assumed  to  accommodate  to  the  gas  properties. 

4-  Cluster- cluster  collisions 

The  cluster-cluster  collisions  is  an  important  process  that  has  to  be  taken  into  account  for  accurate  description 
of  cluster  evolution,  since  it  is  the  key  factor  determining  the  size  distribution  of  larger  clusters.  Cluster- 
cluster  collisions  have  different  outcomes  which  generally  may  be  classified  as  coalescence  and  reflexive 
and  stretching  separations.  The  dynamics  of  water  droplet  collisions  for  macroscopic  particles  was  studied 
experimentally  in  Ref.  38  where  the  boundaries  between  both  of  the  separating  collisions  and  coalescence 
collision  were  examined  as  a  function  of  the  size  ratio  and  the  Weber  number  in  the  wide  range  of  Weber 
numbers  from  1  to  100.  For  microscopic  particles  with  sizes  from  dimers  to  thousand-mers,  the  authors 
are  not  aware  of  any  comparable  to  Ref.  38  systematic  study  where  the  results  of  cluster  collisions  would 
be  analyzed  for  different  Weber  numbers.  Extrapolating  the  results  of  Ref.  38  to  microscopic  particles  of 
interest  in  this  work,  one  can  notice  that  for  typical  plume  temperatures  on  the  order  of  100  K  and  thus 
Weber  numbers  on  the  order  of  the  unity  or  less,  the  clusters  would  mostly  experience  coalescence  and  not 
separation.  However,  such  an  extrapolation,  although  partially  justified  for  hundred-  and  thousand-mers,  is 
much  more  questionable  for  smaller  clusters,  where  more  reflexive  collisions  may  be  expected.  In  this  work, 
in  the  absence  of  reliable  size  and  relative  velocity  dependence  of  the  collision  outcome  for  small  clusters,  a 
constant  coalescence  probability  is  assumed. 

The  cluster-cluster  collisions  are  modeled  using  a  conventional  DSMC  algorithm.  The  majorant  frequency 
scheme^^  of  the  DSMC  method  was  utilized  for  this  purpose.  At  every  time  step,  the  current  maximum 
cluster  size  is  obtained  in  each  cell.  Then,  the  majorant  collision  frequency  is  calculated  based  on  this 
maximum  cluster  size  and  maximum  relative  collision  velocity  evaluated  from  the  local  gas  temperature. 
The  majorant  collision  frequency  is  then  multiplied  by  the  coalescence  probability,  since  only  the  coalescence 
events  are  modeled  (reflexive  separation  is  believed  to  have  negligible  effect  on  cluster  properties,  and  the 
stretching  separation  process  is  not  included  in  the  present  model).  After  a  pair  of  clusters  K  and  L  is 
selected  for  physical  collision,  the  coalescence  event  is  modeled,  with  the  result  being  a  larger  cluster  M 
with  mass  m  and  internal  energy  Eint  calculated  from  the  properties  of  colliding  clusters  using  the  mass  and 
energy  conservation  constraints.  The  laws  dictated 

niM  =  rUK  +  triL,  Eint,M  =  Eint,K  +  Eint,L  —  Q, 

with  Q  =  —Qm  +  Qk  +  Ql,  where  Qi  is  the  energy  of  vaporization  of  cluster  i. 

After  all  collision  and  evaporation  processes  are  simulated  for  a  given  time  step,  the  mass  and  energy 
change  over  this  time  step  are  calculated  over  all  cells  in  order  to  be  included  in  the  Eulerian  gas  flow 
equations.  The  primary  purpose  of  this  step  is  accurate  conservation  of  all  conservative  properties  in  the 
simulation. 
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III.  Thermal  Bath  Relaxation 


Inelastic  cross  sections  for  monomer-monomer  and  monomer-cluster  collision  mechanics  are  needed  as 
part  of  a  comprehensive  validation  of  a  kinetic  condensation  model.  These  cross  sections,  in  general  a 
function  of  the  translational  and  internal  energy  states  of  pre-  and  post-collision  particles,  are  unavailable 
for  the  species  and  temperatures  desired.  However,  equilibrium  rates  for  nucleation  and  evaporation  for  both 
water  and  argon  are  available  in  literature.  Furthermore,  a  necessary  condition  for  the  model  to  fulfill  is 
that  it  produces  correct  behavior  in  equilibrium,  although  it  does  not  guarantee  the  model  will  have  correct 
nonequilibrium  behavior. 

Such  equilibrium  behavior  was  modeled  in  this  work  by  a  thermal  bath  relaxation  of  both  argon  and  water 
at  various  temperature  conditions.  The  equilibrium  constants  for  the  formation  of  clusters  were  calculated 
and  compared  to  published  results  of  Refs.  39,40.  In  addition,  they  were  compared  to  the  previously  obtained 
results  of  the  DSMC  based  model. 

To  model  argon  equilibrium,  over  10  million  simulated  particles  were  used,  and  the  run  was  allowed  to 
run  until  a  steady  equilibrium  value  was  reached,  usually  about  half  a  million  time  steps.  The  timestep 
for  argon  was  selected  to  be  8  x  10“^^  s,  so  that  there  were  much  fewer  than  1  collision  per  molecule  per 
timestep,  making  the  results  independent  of  step  size.  Number  density  of  10^^  m“^  was  selected  to  ensure 
that  clusters  made  up  less  than  0.1%  of  the  gas,  while  maintaining  the  10  million  particle  requirement.  This 
ensured  that  adequate  numbers  of  particles  were  present  for  statistics,  and  that  the  clusters  did  not  have  a 
significant  impact  on  the  behavior  of  the  gas. 

The  dimer  formation  rate  krec  for  argon  was  computed  using  the  stable  dimer  formation  rate  from  Ref.  39, 
that  is  written  as 

krec  =  Ax  exp  —CT  (7) 

The  values  of  A,  B,  and  C  given  in  Ref.  39  are  A  =  10.15  x  10“"^^  m®-molec“^-  s“^,  B  =  —0.278  and 
(7  =  3.10  X  10-3  K-i. 

The  argon  equilibrium  rate  as  a  function  of  temperature  is  shown  in  Fig.  1.  It  is  compared  with  the 
DSMC  results  from  Ref.  30.  There  is  good  agreement  between  this  model  and  the  DSMC  results,  which  were 
shown  to  agree  with  several  experimental  results.  However,  there  is  a  breakdown  in  agreement  at  higher 
temperatures  where  the  present  model  sharply  under-predicts  the  equilibrium  constant.  The  reason  for  good 
agreement  at  lower  temperatures  is  that  the  temperature  dependence  of  the  inelastic  collision  number  Z  was 
chosen  to  allow  this  to  fit  well.  However,  the  under-prediction  at  higher  temperatures  occurs  regardless  of 
the  value  of  Z  chosen.  This  is  acceptable  for  this  work,  as  the  temperatures  of  further  interest  all  fall  within 
the  range  of  good  agreement. 


Figure  1.  Argon  dimer  equilibrium  rate  as  a  function  of  gas  temperature. 

In  modeling  water  equilibrium,  about  1  million  simulated  particles  were  used,  and  the  process  was  again 
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allowed  to  run  until  a  steady  equilibrium  was  reached,  about  half  a  million  timesteps.  The  timestep  was  set 
at  1  X  10“^  s,  which  ensured  that  there  were  still  much  fewer  than  one  collision  per  particle  per  timestep,  but 
also  that  the  system  would  reach  complete  equilibrium  within  a  reasonable  number  of  timesteps.  Similar  to 
argon,  number  density  of  10^^  m“^  was  selected  to  ensure  that  clusters  were  less  than  0.1%  of  the  gas  while 
maintaining  1  million  simulated  particles.  Due  to  the  low  fraction  of  clusters,  this  ensured  the  behavior  of 
the  gas  was  not  influenced  by  the  presence  of  the  clusters.  Note  that  using  a  smaller  gas  density  does  not 
change  the  results  of  the  simulations. 

The  water  equilibrium  rate  is  shown  in  Fig.  2.  It  is  compared  with  the  theoretical  predictions.^*^  There 
is  good  agreement  between  the  current  model  and  the  results  of  Ref.  40  at  all  temperatures  investigated  in 
this  work.  Furthermore,  the  parameter  Z  has  relatively  small  effect  on  the  equilibrium  rate  for  water. 


Figure  2.  Water  dimer  equilibrium  rate  as  a  function  of  gas  temperature. 

The  parameter  Z,  or  inelastic  collision  number  is  effectively  the  inverse  probability  of  the  energy  trans¬ 
fer  between  the  internal  modes  of  a  dimer  and  the  translational  modes  in  the  dimer-monomer  collisions, 
was  found  to  be  an  important  factor  influencing  the  magnitude  of  the  equilibrium  constant.  A  possible 
explanation  is  as  follows.  The  dimers  are  formed  after  three-body  collisions,  and  typically  have  internal 
energies  smaller  than  the  evaporation  energy  after  those  collisions.  In  argon,  the  evaporation  energy  for  a 
dimer  is  small  compared  to  the  typical  total  collision  energy  for  all  the  temperatures  under  consideration. 
{Eevap/k  «  140  K).  This  means  most  dimers  have  an  internal  energy  in  excess  of  the  evaporation  energy 
after  only  a  couple  inelastic  collisions  with  monomers  for  thermal  bath  temperatures  greater  than  140  K. 
With  an  internal  energy  larger  than  the  evaporation  energy,  the  lifetime  of  these  dimers  is  very  short,  about 
a  picosecond.  This  means  the  dimer-monomer  energy  transfer  is  the  main  contributor  to  quick  dimer  disso¬ 
ciation.  It  is  important  to  note  that  Z  has  little  impact  on  dimer  formation,  and  only  affects  the  evaporation 
rates.  Therefore,  in  the  range  of  temperatures  of  interest,  the  equilibrium  constant  for  argon  was  found  to 
be  nearly  proportional  to  Z~^,  except  for  at  high  temperatures. 

For  water,  the  equilibrium  constant  in  much  less  dependent  on  Z .  This  is  because  the  evaporation  energy 
of  a  dimer  is  much  larger  than  the  translational  energy  of  colliding  molecules  and  dimers.  For  example,  the 
reduced  evaporation  energy  for  water  is  Ef.vap/k  «  1,800  K,  while  the  gas  temperatures  were  on  the  order 
of  300  K.  Therefore,  Keq  for  water  is  much  less  dependent  on  Z,  and  is  effectively  independent  of  its  value. 

IV.  Water  cluster  size  distribution  in  nozzle  expansion 

The  second  part  of  the  validation  and  numerical  analysis  of  the  presented  condensation  model  is  focused 
on  the  nucleation  and  evolution  of  small  water  clusters  in  a  conical  nozzle.  The  study  was  prompted  by 
the  availability  of  high  quality  experimental  data**®  on  terminal  size  distribution  of  water  clusters  in  the 
wide  range  of  flow  conditions  where  the  cluster  size  distribution  changes  its  shape  from  exponential  at  low 
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pressures  to  bimodal  at  intermediate  pressure  to  lognormal  at  high  pressures.  The  experimental  results  were 
obtained  by  doping  the  water  clusters  by  one  Na  atom,  which  is  photoionized  close  to  the  threshold  without 
fragmentation.  The  nozzle  is  a  conical  nozzle  with  a  41°  opening  angle,  a  total  length  of  2  mm,  and  a  throat 
diameter  of  50  /im.  Four  different  stagnation  pressures  were  computed,  considered  in  Ref.  18,  1.577  bar, 
2.173  bar,  5.144  bar,  and  8.307  bar,  with  the  corresponding  stagnation  temperature  of  495  K.  Since  the 
background  pressure  effect  in  the  experiment  is  believed  to  be  small, the  flow  expansion  into  the  vacuum 
is  modeled. 

The  computations  were  conducted  on  a  500x150  spatial  grid,  with  cell  sized  reduced  in  radial  direction 
and  increased  in  the  axial  direction.  An  adiabatic  nozzle  nozzle  wall  condition  was  used,  although  previous 
studies  have  found  that  there  is  no  impact  of  the  wall  condition  on  the  coreflow  were  the  cluster  sizes 
are  recorded.  The  number  of  simulated  droplets  was  about  500,000,  which  was  found  to  provide  adequate 
statistical  accuracy  for  the  calculations.  The  particles  were  assumed  to  condense  on  the  nozzle  surface. 
Uniform  inflow  conditions  were  imposed  at  the  nozzle  throat,  calculated  from  the  isentropic  flow  relations. 
To  compare  the  cluster  size  distributions  with  the  terminal  distributions^®  measured  far  downstream  from 
the  nozzle,  the  computed  the  size  distributions  in  several  stations  along  the  nozzle  axis  were  analyzed  to 
provide  distance-independent  distributions.  The  domain  size  was  increased  in  the  axial  direction  from  4  mm 
for  the  lowest  pressure  to  20  mm  for  the  highest  pressure  to  ensure  that  the  size  distributions  at  the  exit 
boundary  are  essentially  frozen. 

Typical  run  time  for  the  lowest  pressure  under  consideration  was  several  hours,  and  for  the  highest 
pressure  was  up  to  two  days  on  a  single  processor  computer.  Comparing  these  numbers  with  those  of  Ref.  30 
where  a  DSMC  method  was  used  to  model  a  1.577  bar  water  expansion,  the  new  approach  is  about  50 
times  faster  than  the  DSMC  based  one  for  the  lowest  pressure,  and  this  factor  will  grow  significantly  with 
pressure.  The  reduction  in  speed  is  mostly  related  to  the  time  efficient  modeling  of  gas  transport  with 
a  continuum  method.  Since  clusters  comprise  only  a  relatively  small  fraction  of  the  particles  in  the  flow, 
gas  transport  modeling  is  the  most  time  consuming  part  of  any  DSMC-based  technique.  Note  that  species 
weights  for  cluster  species  would  reduce  the  time  requirements  of  a  DSMC-based  condensation  model,  but 
the  application  of  weights  is  questionable  in  condensing  flows  since  the  condensation  significantly  changes 
the  gas  flow. 

Consider  first  the  gas  and  particle  properties  along  the  nozzle  axis.  The  gas  translational  temperature 
for  the  lowest  and  highest  pressures  under  consideration  is  shown  in  Fig.  3.  Here,  X=0  corresponds  to  the 
nozzle  throat.  As  expected,  the  water  nucleation  results  in  noticeable  increase  in  gas  temperature  due  to 
evaporation  heat  release.  For  the  1.577  bar  case,  the  temperature  in  the  plume  region  is  up  to  30  K  higher 
when  the  condensation  is  included,  which  is  comparable  to  the  magnitude  of  the  temperature  in  the  non¬ 
condensing  flow.  Note  a  small  temperature  increase  at  about  0.25  mm  is  related  to  the  compression  wave 
that  originates  near  the  nozzle  throat  and  propagate  to  the  nozzle  axis.  It  presents  both  in  the  condensing 
and  non-condensing  flow,  and  the  location  is  nearly  the  same  since  the  impact  of  the  condensation  is  not 
very  significant  at  this  point. 


Figure  3.  Gas  temperature  profile  along  the  nozzle  axis  for  po  =  1.577  bar  (left)  and  po 


8.307  bar  (right). 
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For  the  8.307  bar  case,  the  influence  of  the  condensation  is  obvious  almost  immediately  after  the  nozzle 
throat  (the  temperatures  start  to  deviate  after  the  first  100  /im  from  the  throat),  and  in  the  plume  the  gas 
temperature  is  several  times  higher  in  the  condensing  flow.  Note  some  statistical  scatter  seen  in  this  figure, 
where  instanteneous  gas  properties  are  presented  (an  interpolation  procedure  was  used  here  to  smooth  the 
results).  The  instanteneous  properties  are  dependent  on  current  cluster  properties,  and  the  use  of  a  finite 
number  of  simulated  clusters  contributes  to  their  scatter.  The  higher  temperatures  in  the  nozzle  for  the 
condensing  flow  cause  the  formation  of  a  compression  wave  near  the  nozzle  lip,  that  propagates  downstream 
and  reflects  at  the  axis  at  X«12  mm.  This  results  in  a  significant  rise  in  gas  temperature. 

Consider  now  the  terminal  cluster  size  distributions  at  different  stagnation  pressures.  Note  first  that  there 
are  several  important  properties  that  strongly  affect  the  size  distributions,  among  which  are  the  evaporation 
heat,  heat  capacity,  monomer  sticking  and  cluster  coalescence  probabilities.  The  first  two  of  these  properties 
are  mostly  functions  of  the  cluster  size,  and  the  latter  two,  being  characteristics  of  binary  collisions,  depend 
on  the  cluster  sizes,  internal  energies,  and  relative  collision  velocities.  The  use  of  a  constant  coalescence 
probability  in  this  work  is  a  significant  oversimplification  of  the  actual  cluster  collision  process,  primarily 
related  to  the  lack  of  information  on  collisions  of  small  clusters.  While  the  coalescence  probability  of  two 
relatively  large  clusters  (100-mers  and  larger)  may  be  reasonably  assumed  to  be  close  to  the  unity  for  Weber 
numbers  on  the  order  of  1,  the  coalescence  of  smaller  clusters  is  less  likely  and  for  the  limiting  case  of  dimer 
collisions  may  approach  that  of  monomer  sticking,  which  is  about  0.2  for  water. 

The  numerical  analysis  has  shown  that  the  size  distribution  significantly  depends  on  the  coalescence 
probability,  as  shown  in  Fig.  4.  The  increase  in  the  coalescence  probability  from  0.25  to  1  results  in  significant 
redistribution  of  cluster  sizes  and  shift  from  smaller  sizes  to  larger  ones.  Such  a  trend  is  expected,  since 
higher  coalescence  probability  at  a  given  collision  rate  increases  the  population  of  large  clusters.  Although 
the  coalescence  is  accompanied  by  energy  release  from  the  electronic  structure  of  smaller  clusters  to  the 
internal  energy  of  larger  clusters,  the  larger  internal  energy  is  then  redistributed  over  a  significantly  larger 
number  of  internal  degrees  of  freedom.  The  resulting  gas  temperatures  were  therefore  found  not  to  change 
noticeably  with  the  coalescence  probability. 


Figure  4.  Computed  cluster  size  distributions  for  different  coalescence  probabilities  P  at  po  =  1.577  bar  (left) 
and  po  =  2.173  bar  (right). 

Comparison  of  the  computed  and  experimental  cluster  size  distributions^®  for  these  pressures  is  presented 
in  Fig.  5.  Since  the  dimers  were  not  measured  in  Ref.  18,  hereafter  the  experimental  points  were  normalized 
to  match  the  population  of  the  computed  clusters  excluding  dimers.  The  case  with  a  coalescence  probability 
of  0.5  gives  better  agreement  with  the  data,  and  is  therefore  shown  here.  It  needs  to  be  mentioned  here 
that  a  local  minimum  observed  for  6-mers  and  a  local  maximum  observed  for  8-mers  is  not  a  statistical 
fluctuation,  but  the  consequence  of  the  corresponding  minimum  and  maximum  in  the  cluster  evaporation 
energies.  Note  that  for  po  =  1.577  bar,  the  best  agreement  with  the  data  would  produce  a  computation 
that  utilizes  a  constant  coalescence  probability  between  0.25  and  0.5,  whereas  for  po  =  2.173  bar,  a  larger 
coalescence  probability  between  0.5  and  1.0  would  produce  a  better  agreement.  This  is  reasonable,  since 
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higher  pressures  are  generally  characterized  by  higher  degree  of  nucleation  and  larger  cluster  sizes,  for  which 
the  coalescence  probability  is  expected  to  increase. 


Cluster  Size 


Figure  5.  Terminal  cluster  size  distributions  for  po  =  1.577  bar  (left)and  po  =  2.173  bar  (right):  comparison  with 
data.^® 

For  the  two  largest  pressures  under  consideration,  the  computations  with  a  coalescence  probability  of  1 
provide  better  agreement  with  the  data,  and  the  corresponding  results  are  shown  in  Fig.  6.  For  po  =  5.144  bar, 
the  computed  location  of  the  second  maximum  in  the  distribution  function  agrees  well  with  the  corresponding 
experimental  value,  although  the  population  of  such  clusters  is  somewhat  higher  in  the  experiments.  The 
most  noticeable  difference  is  observed  in  the  large  cluster  tail,  where  clearly  more  clusters  were  observed  in 
the  experiment.  In  the  calculation,  the  large  cluster  tail  is  closer  to  the  lognormal  shape.  Interestingly,  the 
situation  is  opposite  for  po  =  8.307  bar,  for  which  the  tail  is  somewhat  more  populated  in  the  numerical 
prediction.  More  importantly,  the  numerical  results  do  not  produce  a  clear  bimodal  structure  at  this  pressure. 
Although  this  is  clearly  related  to  some  approximations  used  in  the  model,  more  research  is  needed  to  single 
out  the  most  important  reason  for  this.  The  average  cluster  sizes  for  the  above  computation  vs  experiment 
comparisons  are  summarized  in  Table  IV.  There  is  a  reasonable  agreement  between  the  results,  especially 
for  the  three  lowest  pressures. 


Figure  6.  Terminal  cluster  size  distributions  for  po  =  5.144  bar  (left)and  po  =  8.307  bar  (right):  comparison  with 
data.^® 
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Stagnation  pressure,  bar 

Computation 

Experiment 

1.577 

12 

9 

2.173 

18 

20 

5.144 

107 

80 

8.307 

417 

338 

Table  1.  Averate  computed  and  measured  cluster  sizes  at  different  pressures. 


V.  Conclusions 

New  method  for  modeling  homogeneous  condensation  is  presented,  based  on  Eulerian  description  of  gas 
phase  coupled  with  a  Lagrangian  approach  to  the  cluster  phase.  A  continuum,  Euler  /  Navier-Stokes  solver 
VAC  is  used  to  model  the  gas  transport,  and  a  kinetic  particle  solver  is  developed  in  this  work  to  simulate 
cluster  nucleation  and  growth.  The  conservation  of  properties  is  enforced  through  a  two-way  coupling,  with 
gas  properties  influencing  the  cluster  evolution  through  the  dimer  formation  and  monomer-cluster  collisions 
(both  elastic  and  inelastic),  and  mass,  momentum,  and  energy  transfer  from  the  cluster  to  the  gas  phase 
conducted  through  source  terms  in  the  continuum  equations.  The  proposed  approach  is  orders  of  magnitude 
faster  than  a  comparable  approach  based  on  the  DSMC  method.  Note  also  that  it  may  easily  be  extended 
to  model  heterogeneous  condensation. 

The  following  cluster-related  processes  are  taken  into  account  in  the  kinetic  solver:  (i)  collisional  dimer 
formation  that  uses  theoretical  temperature-based  dimer  formation  rates  defining  the  number  of  dimers 
created  in  each  cell  per  time  step,  (ii)  elastic  monomer-cluster  collisions  that  change  the  translational  and 
internal  energies  of  colliding  particles,  with  energy  transfer  modeling  using  the  Larsen-Borgnakke  model, 
(iii)  inelastic  monomer-cluster  collisions  that  result  in  monomer  sticking,  (iv)  cluster-cluster  coalescence 
simulated  with  a  conventional  DSMC  collision  algorithm  based  on  the  majorant  frequency  scheme,  and  (v) 
evaporation  of  monomers  from  clusters  based  on  the  RRK  model. 

The  new  model  was  found  to  reproduce  well  the  known  theoretical  dimer  formation  equilibrium  constants 
for  two  gases  under  consideration,  argon  and  water.  Water  nozzle  expansion  was  modeled  with  the  stagnation 
pressure  ranging  from  1.5  bar  to  8.3  bar,  which  corresponds  to  the  average  cluster  size  increasing  from 
below  10  to  over  300.  The  results  on  the  terminal  cluster  were  found  sensitive  to  the  cluster  coalescence 
probability,  with  the  average  cluster  size  increasing  significantly  when  this  probability  was  increased  from 
0.25  to  1.  Comparison  with  available  experimental  data  have  shown  good  agreement  at  lower  pressures, 
and  somewhat  worse  agreement  at  the  highest  pressure  under  consideration,  where  no  visible  bimodal  size 
distribution  structure  was  noticed  in  the  calculations. 
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